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Free motion around black holes with discs or rings: 
between integrability and chaos — II 



O. Semerak and P. Sukova* 

Institute of Theoretical Physics, Faculty of Mathematics and Physics, Charles University in Prague, Czechia 



(N 

o 

(N 

> 
O 



> 
o 



ABSTRACT 

We continue the study of time-like geodesic dynamics in exact static, axially and reflection 
symmetric space-times describing the fields of a Schwarzschild black hole surrounded by thin 
discs or rings. In the previous paper, the rise (and decline) of geodesic chaos with ring/disc 
mass and position and with test particle energy was revealed on Poincare sections, on time 
series of position or velocity and their power spectra, and on time evolution of the orbital 'lat- 
itudinal action' . In agreement with the KAM theory of nearly integrable dynamical systems 
and with the results observed in similar gravitational systems in the literature, we found orbits 
of very different degrees of chaoticity in the phase space of perturbed fields. Here we compare 
selected orbits in more detail and try to classify them according to the characteristics of the 
corresponding phase-variable time series, mainly according to the shape of the time-series 
power spectra, and also applying two recurrence methods: the method of 'average directional 
vectors', which traces the directions in which the trajectory (recurrently) passes through a cho- 
sen phase-space cell, and the 'recurrence-matrix' method, which consists of statistics over the 
recurrences themselves. All the methods proved simple and powerful, while it is interesting 
to observe how they differ in sensitivity to certain types of behaviour. 
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1 INTRODUCTION 

The black holes discussed in most university courses are isolated, 
stationary and living in asymptotically flat space-times, that is to 
say, they belong to the Kerr(-Newman) family. The metric which 
describes this family is not entirely simple, but it has a number of 
"nice" properties. Among others, its multipole structure — namely 
that required for the black-hole uniqueness theorems to work — is 
just the one that permits the solution of (el ectro-)geo desic equa- 
tions in terms of separated first integrals (e.g. I Will 20091) . This full 
integrability is mostly lost if any of the assumptions (isolation, sta- 
tionarity, asymptotic flatness) is released. More precisely, the inte- 
grability holds in t he whole fam ily of Kerr-Newman-NUT-(anti-)de 
Sitter space-times dCarterlll968h : it is connected with the existence 
of an irreducible second-order Killing tensor (Wal ker & Penrose] 
1 19701) and has been shown to follow from the existence of a princi- 
pal (c onformal) Killing- Yano tensor there (e.g. lFrolov & Kubiznakl 
2008). All these space-times are of Petrov type D and represent sub- 
class of the Plebahski-Demiahski solutions with non-accelerated 
sources. Besides mass, electric charge and rotational angular mo- 
mentum, they also include cosmological constant, magnetic charge 
and NUT parameter (of which the last two do not seem to be phys- 
ically relevant, however). 

Although the Kerr(-Newman) metric is being referred to when 



E-mail: lvicekps@seznam.cz 



speaking about black holes in galactic nuclei and X-ray binaries, 
the above assumptions can only be valid approximately in such as- 
trophysical circumstances; strictly speaking, they are all violated. 
Indeed, the observability of the supposed black holes alone implies 
that they have to be interacting with matter, thus non-isolated and 
non-stationary. Black holes certainly dominate the gravitational po- 
tential and intensity in their wide surroundings, but higher deriva- 
tives of the field (curvature) may be affected by nearby matter sig- 
nificantly. And these higher derivatives govern stability of motion. 
Hence, due to its own gravity, the matter may in fact settle down, 
around a central black hole, to a different configuration than which 
would be assumed by a test (non-gravitating) matter. 

In a non-linear theory like general relativity it is not simple 
to specify what violation of the given (Kerr) metric is already large 
enough to invalidate various related conclusions and to bring physi- 
cal differences with observable consequences. If the ambient matter 
is dilute and its self-gravitation effects are correspondingly weak, 
or when the source is more concentrated but only the field farther 
from it is relevant, the problem is usually being addressed by per- 
turbation techniques. However, perturbation solutions are given in 
terms of series that practically must be truncated somewhere, so 
they cannot represent fully the self-gravitation effects (and in linear 
order they do not encompass them at all). The perturbative descrip- 
tion is mainly questionable in the case of two- or one-dimensional 
additional sources (like discs or rings), because even if the total 
mass of such sources is very small, they however constitute a sin- 
gularity and thus cannot be considered weak in their vicinity. 
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It is sure in any case that almost any deviation from the Kent- 
Newman) metric implied by the presence of additional matter leads 
to the loss of complete (electro-)geodesic integrability. This is even 
true when the additional source keeps the symmetries of the "orig- 
inal", pure black-hole space-time, i.e. when their resulting "super- 
position" is stationary, axially and reflection symmetric and orthog- 
onally transitive (this last property is ensured if and only if the 
source elements do not perform any other motion than steady orbit- 
ing along the direction of axial symmetry). The loss of integrability 
in turn means that motion in the field of such a system is chaotic 
in general. The robustness of this conclusion makes the subtle phe- 
nomenon of chaos one of the aspects of astrophysical black-hole 
systems that should be taken i nto account and that can have o b- 
servable consequences (see e.g. lLukes-Gerakopoulos et al.| [2010). 

Needless to say, in astrophysical systems with accreting black 
holes the matter elements have many other (and more serious) rea- 
sons why to behave in a chaotic way. They extend from micro-scale 
processes over (magneto)hydrodynamics of accretion to interaction 
with (chaotic) radiation. However, as opposed to the case of dis- 
crete sources moving in the field of accreting stellar-mass black 
hole (e.g. in an X-ray binary), these "physical" reasons for chaos 
should be less important in the case of whole stars orbiting a su- 
permassive black hole in a galactic nucleus. Under the presence of 
a heavy accretion disc (and/or massive gas torus farther away), the 
motion of such "test particles" may exhibit chaotic features on a 
sufficient time-scale, due to the "gravitational" reasons alone. Ad- 
mittedly, there are whole star clusters rather than single stars around 
the nuclear black holes, so each individual star would also "feel" 
perturbations from all the other stars. Such a many-body problem 
is very difficult to tackle in general relativity and even in Newto- 
nian theory it is being treated numerically. A possible simplifica- 
tion is to approximate the influence of the whole cluster on one 
particular star by means of an additional spherical or other simple- 
shape potential (see e g. iKaras & Subill2007l ; iMadigan et"ai1l2009l ; 
iLQckmann et alj2009h . 

In spite of the "shadowing theorem' ' dPalmeil 12000), it often 
seems hopeless to model the behaviour of a realistic non-linear dy- 
namical system in detail, because sensitive dependence on initial 
conditions — one of the characteristic signs of chaos — involves sen- 
sitive dependence of conclu sions on our ability to d escribe and treat 
the system exactly (see e.g. ludd & Stemler 2009). Unfortunately, 
the non-linearity of general relativity as the theory of the underly- 
ing configuration space adds another piece of "sensitivity" to the 
problem, and furthermore it strongly limits our compass to treat 
the problem exactly. In particular, the non-linearity severely limits 
the possibility of describing exactly the gravitational field (i.e. the 
configuration space) of multi-component systems, even if they are 
not extended. At present, the exact analytic treatment of black holes 
with additional sources is only practically possible in static and axi- 
ally symmetric case (thus even rotation is excluded in general) with 
zero cosmological constant. Outside of the sources, the complete 
space-time can then be described by the (Weyl) metric containing 
just two unknown functions, one of which has the meaning of New- 
tonian gravitational potential. In a vacuum, the Einstein equations 
yield Laplace equation for this potential (like in the Newtonian de- 
scription), so its "total" form is obtained simply by adding the con- 
tributions from individual sources. The "non-Newtonian" part of 
the problem lies in finding the second unknown metric function by 
a line integral; it is rather an exception than a rule that this can be 
accomplished explicitly. 

We will however not repeat this standard introduction to the 
static and axially symmetric problem; it was summarised (e.g.) in 



the first paper dSemerak & Sukova| [2010). There, we placed un- 
charged annular thin discs without radial pressure (and without heat 
transfer) or their limit — (one-dimensional) rings — symmetrically 
about the (originally Schwarzschild) black hole in order to ap- 
proximate the configuration of matter assumed in black-hole ac- 
cretion systems. Considering, in particular, seve ral discs of the in- 
verted Morgan-Morg an counter- rotating family (Lemos & Letelier 
| l994ISemerakl2003l) an d the Bach- Weyl ring teach & Wevll 19221; 
iD'Afonseca et al.ll2005f) as the additional sources, we studied how 
the dynamics of time-like geodesies in the field of such a system 
depends on parameters, mainly on relative mass and position of the 
external disc/ring and on the energy of test particles. The system 
showed typical features of a weakly non-integrable dynamical sys- 
tem. We observed, on Poincare sections and on time series of phase 
variables and their power spectra, how it gradually turns chaotic 
when relative mass of the disc/ring or energy of the orbits are in- 
creased; we also noticed, quite generically, that for very high values 
of these parameters the system rather recurred to more regular be- 
haviour. 

It is a conventional experience how the originally fully regular 
phase space grains into chains of resonance islands, circumscribed 
by separatrices from which a net of chaotic filaments originates that 
gradually spreads and fuses into a "chaotic sea". However, besides 
the overall development of the phase space, it is also interesting to 
study individual trajectories and try to distinguish different types 
among them. Actually, it is known (KAM theorem) that in "weakly 
perturbed" systems, the phase space contains regions of very differ- 
ent degrees of chaoticity for almost any "strength" of the perturba- 
tion agent. Moreover, even a given single orbit may show different 
degrees of chaoticity/regularity within its different stages. This is 
especially valid for the orbits prone to "sticky motion", i.e. those 
which spend a long time very close to regular islands, while only 
occasionally diverging into a chaot ic sea. We already singled out 
several such orbits in the first paper ( Semerak & Sukova 20 1]3) and 
illustrated that they can be distinguished from "strongly chaotic" 
orbits, drowned in the chaotic sea, according to the power spectra 
of phase-variable time series (that of "vertical" position, z = z(t), 
for example) . In general relativity , this simple method was notably 
employed by Koyama et al. (2007) who demonstrated (on the prob- 
lem of spinning particles in a Schwarzschild background) that the 
power spectra of "strongly chaotic" orbits have "white-noise" low- 
frequency part (relatively flat curve at rather low values), whereas 
the spectra of "weakly chaotic" orbits incline to the "1/frequency" 
shape at low frequencies (and rise to higher values there). This is 
natural since stronger chaos means more irregular time series with 
less distinct periods, in particular without marked recurrences even 
on longer timescales. 

In the present paper, we analyse in more detail several individ- 
ual orbits selected out of those which — at least for a certain time — 
adhere to regular regions ("sticky motion"). We show on Poincare 
diagrams and on power spectra of the corresponding (z-position) 
time series that different parts of these orbits present different de- 
gree of chaoticity. In stages when the orbits tend to fill uniformly 
a large area (chaotic sea), the spectrum approaches white noise at 
low frequencies, whereas almost regular parts of the same orbits 
indeed generate "1/frequency" shape there. However, even "highly 
regular" parts of such world-lines can be clearly distinguished from 
strictly regular orbits. On the other hand, it is known that even 
strongly chaotic deterministic behaviour is clearly distinguishable 
from a random one. It is our object here to check the above for se- 
lected orbits of our dynamical system. We compare how the char- 
acter of geodesic dynamics is revealed by z-position power spectra 
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(section O, by the "average directional vectors" method based on 
statistics over the directions in which the trajectory (rec urrently) 
passe s through specified phase-space cells (section [5J see Sukova 
1201 ll for preliminary results), and by "recurrence plots" which vi- 
sualise the pattern of recurrences to such cells (section^. Several 
useful quantifiers of chaos following from the recurrence plots will 
also be computed and plotted. 



2 POWER SPECTRA OF GEODESICS AND OF THEIR 
PARTS 

Poincare surfaces of section give a picture of how much regu- 
lar/chaotic the system with given parameters is. It is a property 
of "nearly-integrable systems" that some quite irregular trajecto- 
ries already appear after a tiny perturbation, and vice versa, even 
a strongly perturbed phase space still harbours some quite regu- 
lar ones. It can also be recognised on Poincare diagrams that (in 
rather strongly perturbed cases) certain orbits behave quite differ- 
ently within different periods, namely they alternately stick to is- 
lands of regular motion and drift around the "chaotic sea". Hence, 
the "degree of chaoticity" can be ascribed to the system (phase 
space) on the whole, but also to its individual orbits, while in a 
sense it is a property of a given part of a given orbit. It should be 
stressed, however, that such a tracking of regular/chaotic features 
down to a particular segments of particular orbits has to be taken 
with much caution, since it is in fact inconsistent with the essence 
of deterministic chaos as a global phenomenon. This is especially 
clear on systems like billiards where interaction only occurs at dis- 
crete events: one cannot say that their trajectories are almost all 
regular, with chaos solely occurring at this or that point. In systems 
where the interaction is long-range and acts continuously (and, ide- 
ally, even smoothly), like the gravitational one in our case, such at- 
tempts to "localise" the character of dynamics are more propitious, 
but still they do not seem to yield necessary or sufficient criteria 
(cf. lVieira & Letelierll9 96). The global nature of chaos go together 
well with the similar attribute of Fourier transform, which proba- 
bly helps power spectra to be an efficient tool of study of dynamical 
systems. 

We select two geodesic orbits in the field of a black hole (of 
mass M) surrounded by the inverted first Morgan-Morgan disc 
with mass M = 1.3M and inner radius rdi sc = 20M. Such a 
disc is clearly very massive in comparison with what is consid- 
ered astrophysically realistic; it causes quite a strong perturbation 
of the original black-hole field and the geodesic dynamics is rather 
chaotic in general (see paper I). The selected two orbits have spe- 
cific energy and angular momentum at infinity £ = u t — 0.956 
and £ = — 4M, respectively^ and initial conditions very close 
to those of t he orbit drawn in light b lue in figure 16 of the pre- 
vious paper l lSemerak & Suk ova 2010); the power spectrum of its 
z-position is clearly seen in the bottom right panel there, among 
several other "weakly chaotic" orbits showing the "1/frequency" 
spectral shape. Note that the power spectrum is obtained using the 
discrete Fourier transform 
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1 We use geometrised units in which c = G = 1. The particle's proper 
time is denoted by r and its four- velocity by u^. The Schwarzschild-type 
coordinates (t, r, 9, <f>) are employed, with vertical position z = r cos 9. In 
particular, the time t is tied to the time-like Killing symmetry of space-time. 



where N is the total number of samples. The spectrum thus illus- 
trates which frequencies dominate the particle motion perpendicu- 
lar to the plane where the disc or the ring is placed. Both selected 
orbits were followed for a long time (of about 10 7 M of proper time, 
which means for some 10 4 -j- 10 5 orbital periods) and then divided 
in parts in order to check how these sections appear in Poincare di- 
agram and in the power-spectrum plot. Rather than showing the fig- 
ures obtained for all the parts, let us jus t give several examp les. In 
general, one can repeat the conclusion of lKovama et al.ld2007l) . also 
supported in our previous paper: the orbits which appear "strongly 
chaotic" in Poincare diagram (filling rather uniformly thick lay- 
ers there) produce irregular time series whose power spectra have 
white-noise character at low frequencies (rather flat at low level, 
without distinct peaks), whereas the orbits which appear "weakly 
chaotic" in Poincare diagram (remaining close to regular islands for 
considerable periods) produce almost regular time series in the pe- 
riods of "sticky motion" (while irregular elsewhere), which brings 
more power into low frequencies and the spectrum assumes 1// 
shape there. 



In figure Q] the two selected orbits are represented as a whole, 
first on Poincare sections showing passages through the equatorial 
plane in (r, u T ) axes, then on the u r (t) behaviour and finally on the 
power spectra of z(t) evolution. The Poincare sections show that 
the first (left) trajectory spends some time in the vicinity of the split 
primary island, but generally it is rather chaotic. The second (right) 
trajectory fills the chaotic sea less densely and apparently prefers to 
stay in the layers adjacent to the central regular regions. The time 
series z(t), u r (t) of the second trajectory really contain longer pe- 
riods of almost regular oscillations, and also the power spectra of 
z(t) confirm this difference in the above described sense: the sec- 
ond (right) spectrum is a bit less "concave" than the first (left) one, 
it tends more to the power-law, 1// shape, mainly at low frequen- 
cies (where it also ends at somewhat higher values). However, both 
trajectories contain rather chaotic as well as rather regular phases; 
we give two examples for each of the orbits in figures [2] and [3] 
The sections shown on the left are confined to the vicinity of reg- 
ular islands for considerable intervals, whereas those on the right 
are quite chaotic as is clear from the respective Poincare diagrams, 
u r (t) courses and power spectra of z(t) evolutions. 



Several more specific features can be recognised in the figures. 
Comparing the weakly chaotic parts of the orbits, one observes that 
their rather "l//-shaped" power spectra can differ in slope and in 
noise degree, mainly in the high-frequency part. (Sure: it is impor- 
tant, among others, what is the periodicity of the island that the 
orbit adheres to.) Some of these spectra can really be almost fitted 
by a straight line; nevertheless, they typically have distinct peak 
in the middle part and a "valley" to the left of this maximum. On 
the other hand, strongly chaotic parts of the orbits provide "cat- 
back", concave spectral shape which cannot be approximated by a 
straight line; the spectra contain less distinct features, in particu- 
lar, less distinct maximum and valley on its left. However, it would 
be misleading to generalise any tendency seen on certain partic- 
ular series of spectra, because mainly the quasi-regular phases of 
motion can involve different amplitudes and periodicities and thus 
bring different spectral features; this is already clear from the time 
series themselves. In spite of it, we try to put together a series of 
orbital sections with different content of chaos/regularity and illus- 
trate how it proves in power spectra (figure |4). 
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Figure 1. Two geodesic orbits in the field of a black hole (of mass M) surrounded by the inverted first Morgan-Morgan disc with mass Ai = 1.3M and inner 
radius r<jj ac = 20M . Both geodesies have specific energy and specific angular momentum at infinity £ = 0.956 and I = 4Af. The orbit on the left is clearly 
more chaotic, the one on the right spends more time in the vicinity of central regular islands in the Poincare diagram (showing passages across the equatorial 
plane, top row). The time series of u r (given by values recorded at passages through the equatorial plane) are plotted in the middle row and the power spectra 
of the z(t) evolution are plotted in the bottom row. 



3 KAPLAN & GLASS' METHOD OF DIAGNOSING THE 
DEGREE OF CHAOS 

There exists a large number of methods how to recognise the degree 
of chaoticity — or even stochasticity — present in the time evolution 
of a given system. We will not give any comparative analysis, and 



not even a listing of them, but will just illustrate that processing 
of time series in some other way can yield information which is 
interesting to compare with what is revealed by power spectra. As 
an example of a simple but useful m ethod, we will consider the one 
suggested by Kaplan & Glass! dl992l) . It was designed to distinguish 
between deterministic and random systems, but we will see that it 
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Figure 2. Example of a rather regular (left column) and rather chaotic (right column) phases of the orbit shown in the left column of figure [T] Equatorial 
Poincare section (top), u r (t) time series recorded at passages through the equatorial plane (middle) and power spectrum of the z(t) evolution (bottom) are 
shown. 



is quite sensitive and also able to recognise how much chaotic the 
(deterministic) system is. 

The method of Kaplan & Glass is based on monitoring the 
evolution of tangent to the trajectory in small subsets of phase 
space. However, the latter is not the "original" phase space of the 
system (this may be unknown), but its d-dimensional embedding 
"reconstructed" from a given data series x(t) by taking the de- 
layed replicas x(t), x(t — At), x(t — 2Ar), . . ., x(t — dAr) as 



its axes (At is some real time shift ); the reconstru ction is justified 
by the delay embedding theorem of lTakensl Jl98ll) . Such a space is 
then coarse grained into a grid of m d cubes and average directions 
of passages through each of these cubes are added up. Namely, first 
the average direction of each (fc-th) pass of a trajectory through the 
given (j'-th) box is recorded as a unit form v^j of the vector con- 
necting the point where the trajectory entered the box with the point 
where it left it. Then the vectors obtained from a large number (%) 
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Figure 3. Example of a rather regular (left column) and rather chaotic (right column) phases of the orbit shown in the right column of figure \T\ Equatorial 
Poincare section (top), u r (t) time series recorded at passages through the equatorial plane (middle) and power spectrum of the z(t) evolution (bottom) are 
shown. 



of transits through the j-th box are summed (vector addition) and 
the length of the resulting vector Vj is normalised by rij , 



Vi = Wi 



i j 

-E 



on n is checked. For random data, L„ decreases with n roughly as 
in particular, the average displacement per step for random 



-1/2. 



walk in d dimensions is (for large n) given by 



r(f) /T 

r(4\ V nd ' 



(2) 



Finally, the resulting norm Vj is averaged over all boxes which 

were crossed n-times, and the dependence of this average (= L^) where F(.) denotes the gamma function. (Note that 
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Figure 4. The u r (t) behaviour (left column) and power spectra of z(i) evolution (right column) of four different sections of orbits shown in figure[T] taken 
from top to bottom, the 1st, the 2nd and the 4th sections are parts of the 1st orbit (shown in the left part of figure[T}, while the 3rd section belongs to the 2nd 
orbit (right part of figure[T}. The degree of chaoticity clearly grows from top to bottom. 
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Figure 5. A: Power spectra of 6 different parts of the two orbits studied in the last section. In the order from top left to bottom right, their degree of chaoticity 
increases (but here on spectra it is barely obvious for the first 3 or 4 of them). The Kaplan-Glass function A(At) calculated for the same 6 orbital phases is 
shown in part B of this figure. 



n -i/2 •) p or a deterministic system, the average decreases more 
slowly or even remains close to a maximal value of one; this is due 
to the fact that in every small neighbourhood the tangent vectors 
from all transits are almost parallel so the norm of their vector sum 
is almost maximal. (It depends on box size, however: theoretically, 
in the limit of infinitesimally fine grain, which is only conceivable 
for infinitely long data series, L„ = 1 for the deterministic dynam- 
ics.) 



Besides the size of the lattice boxes, the result clearly depends 
on the dimension d and the time lag At. In particular, the choice 
of At may be a subtl e issue, as discussed in the original paper 
(Kaplan & Glass 1992). In order to analyse the dependence of L„ 
on At, one can compare the value of Vj with R%, for each box and 
average this over all occupied boxes, 
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Figure 5 - continued B: The Kaplan-Glass function A(Ar) calculated for the same 6 orbital phases whose power spectra are shown in previous figure[5]A. 
The behaviour of A(Ar) clearly distinguishes between these examples, in particular between the first four of them which appear quite similar according to 
spectra. On the other hand, the last two — quite chaotic — cases are similar here, although their power spectra show somewhat different degree of chaoticity. 



A(At) = 



, l-(^) 2 



(3) 



In a theoretical limit, A = for a random walk, whereas A = 1 
for a deterministic system. In practice, with finite series and finite 
"pixel" size, A falls off roughly as autocorrelation function for ran- 
domised signal, while more slowly for a deterministic signal. 

For our system of geodesic dynamics in the static and axially 



symmetric field of a black hole surrounded by a disc or a ring, we 
choose d = 3 and construct the "phase space" out of the time se- 
ries of the test-particle vertical position z = r cos (r and 9 are 
Schwarzschild radius and latitude), i.e. the space is spanned by the 
axes z(r), z(t-At) and z(r-2Ar). We choose its "elementary- 
grain" size of the order of M (thus volume ~ M s ) and focus on 
the behaviour of A with increasing At. Computation of this de- 
pendence for a large number of orbits (or their parts) confirmed 
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that the method of Kaplan & Glass is able to reveal the degree of 
their chaoticity. For rather regular ("sticky") evolutions, A really 
remains very close to one except when At is just multiple of some 
important orbital period. On the other hand, for strongly chaotic 
(parts of the) orbits, A has a peak almost approaching one for a cer- 
tain value of At (which is smaller than a period and is related to 
the first zero of the autocorrelation function) and than it decreases 
quite rapidly (with some oscillation), because around such orbits 
the geodesic flow is sensitive to initial conditions and the determin- 
istic connection between the position x(t) and x(t — At) fades 
away quickly with the growing time lag. The rate of decrease of 
A(At) is a simpler (more easily comparable) indicator of the de- 
gree of chaoticity of the orbit than the power spectrum. Namely, 
it is just value of A(Ar) what is important, whereas the spectrum 
mainly bears its information in the overall shape and the "degree of 
noisiness", of which especially the latter may be difficult to judge 
and compare. At the same time, the Kaplan-Glass A(Ar) indicator 
seems to be quite reliable, and also "fine" in the sense that it can 
distinguish between evolutions which appear very similar (or just 
not easily comparable) on spectra. 

Let us illustrate this on several parts of our two orbits dis- 
cussed in previous section. Figure [5] shows spectra of 6 different 
orbital sections, of which the 1st and the 5th are parts of the sec- 
ond orbit and the rest are parts of the first orbit. The sections are 
placed (top left, top right, middle left, etc.) in the order of increas- 
ing chaos. Figure [5jB shows the A(Ar) functions for the same 6 
orbital sections. It is seen how the A(Ar) behaviour safely distin- 
guishes between the first four examples which are all "rather regu- 
lar" and have just slightly different character of spectra. (But after 
the clue is provided by A(Ar), one admits that the spectra also 
reveal the same tendency.) On the other hand, A(Ar) appears not 
to be so sensitive in more chaotic regions: starting from a certain 
amount of chaoticity, it does not decrease any more and gives sim- 
ilar behaviour for orbital phases which still quite differ in spectra, 
mainly in the low-frequency tail. (For instance, the last two exam- 
ples of figure |4]yield the same course of A(Ar).) 

We add two more illustrations, one (figure |6]l showing two or- 
bits in the field of a black hole with the inverted Morgan-Morgan 
disc of different parameters than above (mass M = 0.5A/ and 
inner radius rdi sc = 14M), and the other (figure [7J showing two 
orbits in the field of a black hole surrounded by a thin, Bach-Weyl 
ring (of mass M = 0.5M and radius r r i ng = 2QM ). The orbits in 
figure [6] are only slightly different, one belonging to five-periodic 
regular islands and the other sticking closely to these islands; you 
can see how they differ in power spectra and in the A(At) depen- 
dence. The orbits presented in figure|7]differ much more, one lying 
very close to the central circular orbit of the primary regular island, 
while the other filling the chaotic sea around and in the vicinity of 
the ring. 



4 RECURRENCE ANALYSIS OF THE ORBITS 

Let us stress again that the above method is just an example of a 
simple but apparently quite reliable way of judging and classify- 
ing the regularity / chaoticity / stochasticity of a given time series. 
Another example of a simple and powerful method is the recur- 
rence analysis which is based on checking the recurrence of or- 
bits of a dynamical syste m to a chosen (small) cell of phase space 
(see lMarwan et alj|2007l for a thorough survey). The pattern of re- 
currences encodes quite credibly the character of the dynamics, 
clearly distinguishing between regular, chaotic and random evolu- 



tions. Within g eneral relativity, the m ethod has recently been em- 
ployed e.g. by Kopa cek et ail ( bold) in their study of orbital dy- 
namics of a charged test particles around a rotating black hole im- 
mersed in a magnetic field. 

A very useful tool for vi sualisation of the recu rrences are re- 
currence plots, introduced bv lEckmann et al.l dl987l) . Suppose one 
knows the phase-space trajectory with constant step of time (we use 
proper time in our case of geodesic motion in a given space- time)0 
Denoting by Xi — X(ji) the N successive points of the phase 
trajectory, one defines the so called recurrence matrix by 

^(e) = e(e-||X i -X J ||), i,j = l,...,N, (4) 

where e is the radius of a chosen neighbourhood (it is called thresh- 
old), || ■ || denotes the chosen norm (in accord with a common ex- 
perience, the picture of long-term dynamics only slightly depends 
on which norm is chosen) and O is the Heaviside step function. 
The matrix thus contains only units and zeros and can be easily 
visualised by representing l's by black dots (while 0's by blank 
spaces) at the respective coordinates i, j. For regular systems, the 
black points tend to arrange in distinct structures, in particular in 
long parallel diagonal lines (their distance scales with period) and 
checkerboard structures, whereas for random behaviours the black 
points are scattered without order. Chaotic systems yield the most 
"artistic" plots: they contain blocks of almost-diagonal patterns as 
well as irregular ones, apparently placed one over another within 
horizontal and vertical structures. The main diagonal Ri t i is triv- 
ial ("line of identity") and present in every system (for it is of- 
ten omitted) and the matrix is symmetric with respect to it. The 
almost-regular blocks correspond to time intervals when the trajec- 
tory sticks to some unstable periodic orbit; the more unstable this 
orbit is, the earlier the trajectory deviates from it and so the smaller 
is the block. Horizontal/vertical lines indicate periods when the sys- 
tem is trapped in some region of phase space without much change. 

Judging the prominence of diagonal or other patterns within 
the recurrence plot by pure observation is of course subjective, so 
several quantifiers of the recurrence-matrix properties have been 
proposed. The simplest of them is the recurrence rate, given by 
ratio of the recurrence points (black ones) within all points of the 
matrix, 

JV 
i,i=l 

Another, most important quantifier is the histogram of diagonal 
lines of a certain prescribed length I, 

N 

P(e,l) = ^2[l-R i - 1J - 1 (e)][l-R i+lJ+l ^)]^ 

l-l 

x]Ji?i+fc,j+fc(e). 

k=0 

Several further quantities can in turn be computed from this his- 
togram. The one called DET is given by ratio of the points which 
form a diagonal line longer than a certain l m i n within all the recur- 
rence points, 



2 The knowledge of just one phase variable (e.g. position) is enough in fact, 
since the phase space can be reconstructed from a sequence of its time- 
delayed series as in the directional-vectors method discussed in previous 
section. 



© RAS, MNRAS 000,[j]-?? 



Chaos around black holes with discs or rings 1 1 




♦ ; + + + t * + + * + + + _ 



A 



i + a- + + + t ~ 

+ + + + + + + + ++ + + + ^ + + 

♦ + + x + + ,+ 



+ + 

+ + 



Figure 6. Two geodesic orbits in the field of a black hole (of mass M) surrounded by the inverted first Morgan-Morgan disc with mass M = 0.5M and inner 
radius t^ sc = 1AM . Both geodesies have specific energy and specific angular momentum at infinity E = 0.955 and t = 3.75A/. The first (left column) 
belongs to five-periodic regular islands and the other (right column) sticks to the first rather closely on Poincare diagrams (top row). The slight difference 
between the orbits also shows on power spectra of their z(t) evolution (middle row) and on the behaviour of their tangent's autocorrelation represented by the 
A(A-r) dependence (bottom row). 



DET(e) 



EL^M) 



L(e) = 



EL '^>0 



P(e,l) 



while the average length of diagonal lines is 



The length of the longest diagonal L max (e) = maxi-i..jv{/i} and 
its inverse DIV(e) = 1/L max (e) are also of interest, since they 
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20M. 

Both geodesies have specific energy and specific angular momentum at infinity € = 0.93 and I = 3.75M . The first (left column) lies deep in the primary 
regular island, while the other (right column) fills the chaotic sea around and in the vicinity of the ring (see Poincare diagrams in top row). The big difference 
between the orbits also shows on power spectra of their z(t) evolution (middle row) and on the behaviour of the Kaplan-Glass averaged autocorrelation 
parameter A (At) (bottom row). 



are most directly related to the rate of divergence of nearby orbits 
and serve as a rough estimate of the largest Lyapunov exponent. 

It is possible to make this estimate more precise, because the 
cumulative histogram of diagonals turned out to yield the second- 
order Renyi's entropy K2 (also called correlation entropy) which 
stands for a lower estimate of the sum of positive Lyapunov ex- 



ponents. Let the phase space be divided in boxes with the size e, 
numbered in some order by ii, . . . Denote by Pi lt ... t i l (e) the 
probability that the point X(At) is in the ii-th box, the follow- 
ing point X(2At) is in the 12-th box, and so on, up to the point 
X(IAt). The second order Renyi's entropy K2 is given by 
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1 X — ^ 2 

K2UJ) — — lim lim lim — — In > p;, i, (e) 



(5) 



Realising that the occurrence of (non-trivial) diagonal of length / 
means that I successive points of the trajectory X(t), X(t + At), 
. . ., X(t + (I — 1) At) lie in the e-neighbourhoods of certain other 
I successive points X (tq), X(tq + At), . . ., X (to + (I — l)Ar), 
iMarwan et alj d2007l) showed that (under the assumption of ergod- 
icity) the sum can be approximated by 



(6) 



where pt(e, I) is the probability of finding the line of length / in the 
boxes centered at points X(t), X(t + At), . . ., X(t + (1-1)At). 
Hence, K2 can be estimated from the relation 



K 2 (e,l) ^ K 2 (e,l) = -—ln Pc (e,l) , 



(7) 



where p c (e, I) is the probability of finding a diagonal line whose 
length is at least I. This means that Ki is determined by a slope 
of the cumulative histogram plotted (in logarithmic scale) against 
the diagonal length I. At the same time, the correlation entropy K% 
was shown to yield a lower estimate of the sum of positive Lya- 
punov exponents, so the estimate K2 is a good indicator of chaotic 
behaviour. 

Similarly as with diagonal lines, one can plot the histogram of 
vertical (or horizontal) lines (against their length v), 



P(e,v) 



N 

E 



and also the respective measure of vertical structures (parallel of 
DET, called laminarity) 



LAM(e) 



Y N vPlv) 



The average length of vertical lines, 



TT(e) 



r" p(v) 



is called trapping time, because it indicates for how long the system 
is "trapped" (without evolution in phase-space). Finally, from the 
probabilities that some chosen diagonal or vertical line has length 
I, p(l) = P(l)/Ni and p(v) = P(v)/N v , where iV ; and N v are 
total numbers of diagonal/vertical lines, one can compute the so 
called Shannon entropies 

N 

L ENTROPY = - p(0 ln P(0> 

i=imin 

N 

V ENTROPY = - V p(v)lnp(v). 

v ~ v min 

Also relevant as indicator is the size of the "white gaps" be- 
tween vertical (or horizontal) lines, because it is related to the re- 
currence times. For example, let us choose a point Xi on some 
trajectory and record the sequence of points which fall in its e- 
neighbourhood, {Xj 1 , Xj 2 , . . .} (this set corresponds to black dots 
in a certain column/row of the recurrence matrix). Compute the 
recurrence times given by differences between serial numbers of 



the consecutive recurrence points Xj k+1 , Xj k multiplied by the re- 
spective proper-time steps, Tfc(e) = (jk+i — jk)Ar. The mean of 
Tfe is called the recurrence time of the first type, Tl. However, the 
recorded set usually also contains successive points between which 
the orbit does not leave the given e-neighbourhood, so Tj.(e) = 1. 
These points (called sojourn points) does not represent true recur- 
rences, so they should be discarded from statistics. After removing 
all the sojourn points, the above recurrence set contains just be- 
ginnings of the vertical black lines. An average of their distances 
(average length of white vertical gaps) is called the 2nd-type recur- 
rence time, T2. The recurrence times typically behave inversely to 
RR. 

For the recurrence analysis to yield plausible results, it is cru- 
cial to set the parameters e, At, lmi n , v m - m properly. The depen- 
dence of the results on these parameters is itself interesting to ex- 
plore. In particular, the recurrence matrix and all the quantities 
computed from it depend critically on the "target" size e. For exam- 
ple, when e is chosen too large, the time step At rather small and/or 
£min (or D m in) too small, there are plenty of sojourn points in the re- 
currence matrix. On the other hand, if e is too small, the matrix may 
come out too sparse. The Z m i n limits the occurrence of short diag- 
onal lines which are often formed, when e is large enough, due to 
the fact that the e-neighbourhoods of the n-th and (n + l)-st recur- 
rence loops have finitely long non-empty intersection, even though 
the loops may be quite diverging from each other. Finally, there is 
a subtle issue of misleading long se condary diagonals which was 
pointed out bv lMarwan et alj d2007l) ; let us only touch on it here: 
an abundant occurrence of sojourn vertical sequences involving the 
main diagonal leads to the occurrence of another long diagonals in 
the main-diagonal vicinity, which deform the recurrence quantifiers 
if taken into account. Namely, the usage of w m i n only excludes the 
short sojourn verticals from the statistics, but does not influence the 
statistics of diagonals; and the usage of Z m i n only excludes short 
diagonals, not the long ones. Therefore, in order to get rid of the 
latter, a certain lower bound for the length of vertical lines (called 
Theiler's parameter, w) is sometimes required already in recording 
the recurrence matrix (in addition to the careful choice of e), i.e. 
only the points are considered whose serial indices i and j satisfy 
\j — i\ > w. We used the Theiler's parameter only in computing 
figures[lO]and[T4](right column). 

Several examples of recurrence plots obtained for our black- 
hole-disc system are attached. For computation of the recurrence 
matrix, we used the Euclidean norm (both position and velocity are 
processed as vectors in Euclidean space) and "Cartesian" spatial 
mesh corresponding to spherical Schwarzschild coordinates; each 
of the three coordinates has been normalised to zero mean value 
and unit standard deviation. First, we drew the recurrence plots for 
orbits (or their parts) studied in previous sections in order to com- 
pare the information they reveal with that provided by Poincare 
maps, time series and their spectra, and by the average-directional- 
vector method. Figure [8] shows recurrence plots of the two orbital 
sections (one rather regular and one rather chaotic) compared in fig- 
ure [3] The almost regular section occupies the upper left half and 
the rather chaotic section occupies the lower right half; the contrast 
in their "entropy" is obvious. Figure [9] brings recurrence plots of 
the six orbital sections whose z = z(t) power spectra and Kaplan- 
Glass averaged autocorrelation parameter A(At) were presented 
in figures [5]A and|5]B, respectively. The plots are again ordered 
in the direction of increasing chaoticity. The first of them really 
looks like one of the dark boxes of the subsequent plots (these 
boxes correspond to "sticky-motion" periods when the orbit is tied 
to a regular region); the second plot shows a clear checkerboard 
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Figure 8. Recurrence plots for the two orbital sections whose Poincare maps, z(t) evolutions and their power spectra were shown in figure [3] We take 
advantage of the recurrence-matrix symmetry and give just halves of the recurrence plots in one square: the almost regular section (the left column of figure[3) 
is above the main diagonal (upper left triangle), while the rather chaotic section (the right column of figure|3j is below the main diagonal (lower right triangle). 
The axis values indicate proper time in units of M. 



pattern (the orbital section is rather regular) which then gradually 
give way to horizontal/vertical pattern at times when the orbit's be- 
haviour changes to chaotic|f|Figure[TO]brings recurrence plots and 
diagonal-length histograms for the two rather regular orbits whose 



3 Note that the middle-row plots of this figure are counterparts of the z(t) 
time series and spectra shown in the first two rows of figure [4] One can 
check by comparison that the transition to chaos in the z(t) evolution really 
occurs at times when the recurrence matrix changes from checkerboard to 
horizontal/vertical regime ("carpet edge" in the recurrence plots). 



Poincare diagrams, z(t) power spectra and tangent's autocorrela- 
tion parameter were shown in figure [6] The recurrence analysis 
well distinguishes between the orbits, in particular, the histogram 
slope (indicated in red in the plots) which yields Renyi's entropy is 
steeper for the less regular orbit. 

Figures QT] and [12] show six (the latter only four) of 
the recurrence-analysis quantifiers (RR, DET, DIV, LAM, 
V ENTROPY and T2), computed for large collections of orbits 
ejected from certain radial ranges within the equatorial plane of 
the black-hole-disc systems. The quantifiers evidently distinguish 
between regular and chaotic orbits, but their sensitivity to differ- 
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Figure 9. Recurrence plots for the orbital sections whose z(t) power spectra and Kaplan-Glass parameter A(Ar) were shown in figures[5](A and B). The z{t) 
evolutions and power spectra corresponding to the middle-row plots are given in the first two rows of figure[4] Proper time goes along the axes in units of M. 
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Figure 10. Recurrence plots and histograms of diagonal-line lengths for the two orbits whose Poincare diagrams, z(t) power spectra and tangent's autocorre- 
lation parameter were shown in figure|6] Both orbits are rather regular: left orbit belongs to five-periodic regular islands and the right one sticks to them. Both 
have been followed up to 2.1 ■ 10 6 M of proper time, with 15M step; we set e = 1.1 and employed the Theiler's parameter ui = 4. The recurrence analysis 
clearly distinguishes between the orbits, in particular, the Renyi's entropy K2, read off from the slope of the cumulative histogram plotted (in logarithmic 
scale) against the diagonal length I (bottom plots), comes out about — 1.5 • 10~ 6 for the regular orbit while about —8.27 ■ 10~ 5 for the "sticky" orbit (notice 
that the x-axis ranges are different, so the difference in slopes is much bigger than how it appears at first sight). Obviously the slope has to be determined from 
the middle part of the histogram which really reflects recurrence properties and where it can be approximated by a straight line. 



ent phase-space features is somewhat different. In the figure [TT] 
the difference between regular and chaotic region is more distinct, 
because the orbits in the left part belong to a large regular island 
which then quite immediately passes into a chaotic sea (right part 
of the plots). On the other hand, the figure Q2] scans a more com- 
plex phase-space region where several regular islands and chaotic 
layers are crossed. For a better picture of which region have been 
considered, the respective Poincare diagram is also shown in figure 
1121 In the chaotic regime, the recurrence rate and the fraction of di- 



agonal lines are seen to be lower, whereas the vertical measures T2 
and V ENTROPY grow significantly. This means that the system 
takes more time to get back to the chosen e-neighbourhood and that 
the vertical-line length oscillates more than in the regular regime; 
namely, in the regular regime, there are either no vertical lines at 
all, or they are of approximately the same length (note, for exam- 
ple, that if v m i n is chosen too small, one gets a lot of "sojourn" 
vertical lines of the same length). The DIV parameter of course 
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Figure 11. Six of the recurrence-analysis quantifiers, computed for 400 geodesies ejected tangentially (with u r = 0) from radii between r = 21. 5M and 
r = 23.5M (with step 0.005M) from the equatorial plane of the system of a black hole (of mass M) and the inverted 1st Morgan-Morgan disc with mass 
1.3M and inner Schwarzschild radius 20M . All the orbits have specific energy E = 0.934 and specific angular momentum I = AM. The orbits have been 
followed for about 250000A/ of proper time (thus the minimal achievable DIV is 4 ■ 10 -6 ) with "sampling period" At = 45M, the minimal length of 
diagonal/vertical lines has been set at 90M and the recurrence threshold at e = 1.1. The left part of the plots scans trajectories belonging to a large secondary 
regular island, while the right part goes through the chaotic sea. All the quantifiers clearly distinguish between the two regimes. 
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Figure 12. The quantifiers RR, DET, DIV and V ENTROPY computed, like in figure[TT] for 470 geodesies ejected tangentially (with u r = 0) from 
radii between r = 5M and r = 24M (by 0.04M) from the equatorial plane of the system of a black hole (M) and the inverted 1st Morgan-Morgan disc 
with mass M = 0.5M and inner radius r^ sc = 18M. All the orbits have specific energy £ = 0.9532 and specific angular momentum £ = 3.75M. The 
orbits have again been followed for about 250000M of proper time with "sampling period" At = 45M, the minimal length of diagonal/vertical lines has 
been set at 90Af and the recurrence threshold at e = 1.25. Poincare diagram of the system is shown at the top, with orbits coloured according to the value of 
their DIV (top left) and according to the slope of the diagonal-line histogram (top right); the scale is logarithmic. A system with a more complicated phase 
portrait has been chosen here in order to compare the results with the simple case of figure fTTI and in order to check how sensitive the quantifiers are to more 
tiny phase-space features. 
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Figure 13. Dependence of four recurrence quantifiers on the choice of the threshold approach e. We go through the geodesies starting with u r = from the 
radii r = 5M 4- 45 M from the equatorial plane of the black-hole-disc system with the disc mass M = 1.3M and inner disc radius r^s c = 20M. All the 
geodesies have specific energy E = 0.956 and specific angular momentum I = AM. Their DIV, LAM, V ENTROPY and T2 parameters are computed 
for 1 1 different values of e, namely 0.50, 0.65, 0.80, 0.95, . . . 1.85, 2.00; in this order, the line colour shifts from dark violet, violet, light blue, dark blue, 
to red and dark red. A primary regular island and two smaller islands aside are evident. 



increases in the chaotic region as well as the "laminarity". All the 
quantifiers get considerably more noisy in the chaotic regions. 

In order to illustrate how Poincare diagrams can be supple- 
mented by information provided by recurrence analysis, the pas- 
sage points are coloured according to the values of selected re- 
currence quantifiers which apply to the respective orbits (so the 
colours are of course mixed in chaotic regions); namely, the top 
left diagram is coloured according to the value of DIV and the top 
right diagram according to the slope of the diagonal-line histogram. 
We have chosen exactly these two quantifiers, because DIV is 
one of the simplest, whereas the diagonal-histogram slope is more 
"sophisticated", tightly connected with Lyapunov exponents. How- 
ever, it is seen that the two colourings provide almost the same 
information (only the colour scales are somewhat different, natu- 
rally). Therefore, the DIV quantifier seems to be more suitable for 
a quick distinction between regularity and chaos, mainly if large 
sets of trajectories are to be processed routinely. The cumulative- 
histogram slope is theoretically more sophisticated, but also harder 
to get (it is "higher-level") and, mainly, its "automatic" computer 
evaluation is much more tricky, because it has to be determined 



from a proper part of the histogram. Namely, only a certain middle 
part of the histogram is relevant, since its short-length end typically 
"diverges" due to increasing number of sojourn points, while the 
long-length end typically falls off quickly due to the finite length 
of the trajectory. (The histogram has to be computed in the limit 
I — > co, so the short-length end has actually no sense; the long- 
length end is of course determined by the fact that practically the 
trajectories cannot be infinitely long.) 

Finally, the dependence of several quantifiers on the choice of 
e is illustrated in figure[T3] There, the behaviour of DIV , LAM, 
V ENTROPY and T2 over geodesies starting from a wide range 
of radii is plotted for 1 1 different values of e. In all the plots, one 
can clearly recognise a large primary regular island and two smaller 
ones. The quantifiers generally show monotonous dependence on 
e — LAM and V ENTROPY increase whereas DIV and T2 
decrease with e, as expected. Figure[T3lmav also serve as a further 
support for the DIV quantifier (see figure[T2]as well). As also con- 
firmed by other similar figures we are not showing here, the large 
regular islands can be recognised easily by any of the quantifiers, 
because the latter fluctuate there much less and around markedly 
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Figure 14. Dependence of the diagonal-line cumulative length-histogram (top row) and of the resulting value of the correlation entropy K'z (bottom row) on 
the choice of the threshold approach e, plotted for a regular (left column) and weakly chaotic (right column) geodesies. Both geodesies have been followed for 
about 2 ■ 10 6 M of proper time. Both have specific energy £ = 0.956 and specific angular momentum £ = AM and live in a space-time of a black-hole (M) 
surrounded by the inverted 1st Morgan-Morgan disc with mass A4 = X.3M and inner radius r,ji S c = 20 M. The recurrence matrix of the chaotic geodesic 
has been computed using the Theiler parameter growing linearly with e (from w = 2 to w = 6). The K2 entropy of this orbit is considerably more sensitive to 
e than that of the regular orbit on the left (the range along both ^2 -axes corresponds to a change by a factor of 4), however the two cases (regular/chaotic) are 
clearly distinguished as seen on the orders at the vertical axes. 



different values than in the chaotic regions. But only the DIV 
parameter appears to be fairly indicative of small islands as well. 
Namely, in regular regions (both large and small) it is typically 2 or- 
ders of magnitude lower than in chaotic regions (DIV ~ 3 x 1CP 6 
within regular regions, about an order higher in thin chaotic layers 
and about 3 x 10~ 4 in large chaotic regions; this chaotic-sea value 
corresponds to a divergence time of several thousands M, which 
represents some ten cycles about the black hole). For longer trajec- 
tories (than those we have treated here) the difference in DIV be- 
tween regular and chaotic regions tends to be even bigger. The other 
quantifiers only respond to small islands by reducing their oscilla- 
tion, but not by a noticeable change of the mean value. Loosely 
speaking, their sensitivity profile is shifted towards larger regu- 
lar regions (see the RR quantifier in figure [T2] for example). The 
last figure [14] illustrates the difference between regular and chaotic 
geodesic on a markedly different value of K2 entropy inferred from 
the cumulative length-histogram of diagonal lines, and also on dif- 
ferent dependence on e of this histogram. 



5 CONCLUDING REMARKS 

In paper I, we performed an overall check of the time-like geodesic 
dynamics in the field of a Schwarzschild black hole surrounded by 
an axially symmetric static thin disc or ring, and of the dependence 
of its chaoticity on parameters characterising the additional source 
and the test particles. In the present paper, we have focused on in- 
dividual orbits and their different parts and showed how the de- 
gree of their irregularity, already visible on Poincare diagrams, can 
be judged in more detail by studying the time series obtained for 
phase variables. We have mainly considered the times series z(t) 
of the position perpendicular to the disc/ring plane, computed their 
power spectra and compared the information thus gained with the 
one provided by the method of Kaplan & Glass which tracks auto- 
correlation between different parts of the series in dependence on 
time shift, and also with outcomes of the recurrence-matrix analy- 
sis of Eckmann et al. All these methods prove simple and powerful, 
while they differ in sensitivity to specific types of behaviour. 
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Let us also include a brief ment ion of several others' results 
which appeared in the field recently. Brink d2008h contributed to 
the topic of complete geodesic integrability, mentioned in the intro- 
duction, by studying geodesies in stationary axisymmetric vacuum 
space-times and the correlatio n of its fab r ic wit h the existence of 
the "fourth integral". (Cf. also iMarkakisI J2012t) for a Newtonian 
treatm ent of the integrability problem.) Verhaaren & Hirschmannl 
(2010) returned to the study of dynamics of test particles with 
spin in a Schwarzschild space-time and argued, on Poincare dia- 
grams and Lyapunov exponents, that smaller values of spin than 
previously t hought can already make the particle motion chaotic. 
iKovacs et al.l fcOllh performed the first post-Newtonian analysis 
of the Sitnikov system (motion of a test body in the field of an 
equal-mass binary, along a line perpendicular to the orbital plane 
and going through the barycentre) and demonstrated, by numerical 
study of the system in dependence on gravitational radius of the 
"primaries", that the relativistic effect of pericentr e adva nce does 
not destroy its chaotic aspects. iRamos-Caro et aL studied 
the motion of test particles in the field of a centre with quadrupole 
deformation surrounded by finite thin discs obtained by super- 
positions of members of a counter-rotating Morgan-Morgan fam- 
ily. They found there is a close connection between linear stabil- 
ity/instability of equatorial circular orbits and regularity/chaoticity 
of general three- dimensional orbits p assing through their radii □ 
The same group l lLetelier et al.ll20l"lD also revisited geodesic dy- 
namics in the system composed of a monopole or an isotropic har- 
monic oscillator and obl ate quadrupole , and f ound several new fea- 
tures not noticed before. IWang& WvJteOllh employed a pseudo- 
Newtonian potential in order to superpose a rotating black hole 
with a quadrupole halo. They analysed emission of gravitational 
waves from particles orbiting in such a field and demonstrated that 
the radiated amplitude and power are close ly related to the degree 
of chaoticity of the orbit. Gal avij d20 1 lh considered the evolu- 
tion of a compact binary perturbed by a third body. Within a post- 
Newtonian version of the Hamiltonian ADM formulation, and us- 
ing basin-boundary analysis and Lyapunov exponents, he examined 
the relative importance of d ifferent PN orders i n indu cing chaos 
in the system. Very recently. ICon topoulos et al. (2012) have anal- 
ysed the classical system of two coupled oscillators and free mo- 
tion in the general relativistic Manko-Novikov space-time (which 
describes a rotating axisymmetric compact body); they mainly fo- 
cused on periodic orbits of the syste ms and on depend ence of their 
properties on orbital energy. Finally, llgata et alj J20T 1) observed on 
Poincare maps that the time-like geodesies bound in the field of the 
Emparan-Real 5D black ring show chaotic features. 

To conclude, it should be admitted that once the evolution of 
a given dynamical system is mastered with sufficient numerical ac- 
curacy, it is rather easy to produce various decorative figures. But 
it is also true that these can indeed yield a good picture of how 
much irregular the system is and how the irregularity depends on 
system parameters. This in turn indicates how much the "perturba- 
tions" acting on real systems degrade the corresponding simplified 
exact models and helps to evaluate the validity of various approxi- 
mations. Turning to our particular problem of orbiting in a static 
compact-centre space-time, it would be suitable to employ such 
methods in order to estimate and compare the significance of vari- 
ous "perturbations" present in real astrophysical situations. Besides 



4 We would like to remember professor Patricio Letelier who was a leading 
expert in the fields of general relativity and chaotic dynamics. He left us just 



at the time when the above paper appeared in MNRAS. 



the gravitational influence of additional matter, discussed (in a sim- 
ple, static and axially symmetric case) in the present work, there 
would also occu r mechanical inter action of the orbiter with that 
matter (see e.g. ISubr"& Karas 2005); both the compact centre and 
the matter around would probably have non-zero angular momen- 
tum, so dragging effects should be incorporated; actually the orbiter 
may itself be endowed with spin or even higher moments; incoming 
gravitational waves can perturb the system as well as those emitted 
by the orbiter (back reaction); of course, if the orbiter was charged, 
it would also be affected by electromagnetic field, if there is some 
around. 

When thinking about possible perturbations of the origi- 
nally completely integrable problem of free test motion in a 
Schwarzschild or Kerr field, one has mainly in mind the motion of 
individual stars around supermassive black holes in galactic nuclei. 
As already stressed in Introduction, there are in fact whole clusters 
of stars in galactic nuclei, so when solving the motion of an individ- 
ual satellite, one should also take into account gravity of the whole 
cluster, or solve that motion right as a part of the problem of N in- 
teracting bodies. Such a problem is difficult within general relativ- 
ity, mainly if the black-hole centre and the disc or ring/toroid should 
also be taken into ac count, but it is being considered within Newto- 
nian th eory (see e.g. lSubr et al .l2004l ; lHaasetal.l2011al ;l Haas et al.l 
1201 lbl and references therein) as well as in post - Newtonian and 
post-Minkowskian approximatio n (e.g. IChull2009l ; iLedvinka et al.l 
l2008l;lHartung & SteinhoffeOllh . 

There are several immediate options for further work. One 
can of course check the results with yet other methods (Lyapunov- 
type coefficients, various other "entropies", basin-boundary anal- 
ysis, Melnikov integral, etc.), find and study particular significant 
orbits of the system (periodic and homoclinic/heteroclinic orbits) in 
detail, or compare the relativistic analysis with Newtonian, pseudo- 
Newtonian or post-Newtonian one. However, we would mainly 
like to focus on astrophysically more realistic situations (parame- 
ter ranges) in the following. This should involve, among others, the 
question of whether to incorporate also another gravitating compo- 
nents like spheroidal halo, disc plus outer ring (or thick toroid), 
or/and jets, and — most importantly — the question of how to ac- 
count adequately for rotation. 
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